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We numerically solve the Klein-Gordon equation at second order in cosmological perturbation 
theory in closed form for a single scalar field, describing the method employed in detail. We use 
the slow-roll version of the second order source term and argue that our method is extendable to 
the full equation. We consider two standard single field models and find that the results agree 
with previous calculations using analytic methods, where comparison is possible. Our procedure 
allows the evolution of second order perturbations in general and the calculation of the non-linearity 
parameter /nl to be examined in cases where there is no analytical solution available. 
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I. INTRODUCTION 



Cosmological perturbation theory is an essential tool for the analysis of cosmological models, in particular as the 
amount of observational data continues to increase. With the recent launch of the planck satellite, the WMAP mission 
reaching its eighth year, and a host of other new experiments, we will have access to more information about the early 
universe than ever before [TJ |5] . 

To distinguish between theoretical models it is necessary to go beyond the standard statistical analyses that have 
been so successful in the recent past. As a result much interest has been focused on non-gaussianity as a new tool to 
help classify and test models of the early universe. Perturbation theory beyond first order will be required to make the 
best possible use of the data. In this paper we outline an important step in the understanding of perturbation theory 
beyond first order, demonstrating that second order perturbations are readily amenable to numerical calculation, even 
on small and intermediate scales inside the horizon. 

Inflationary model building has for the past few years focused on meeting the requirements of first order perturbation 
theory, namely that the power spectra of scalar and tensor perturbations should match that observed in the Cosmic 
Microwave Background (CMB). Inflationary models are classified and tested based on their predictions for the power 
spectrum of curvature perturbations, the spectral index of these perturbations and the ratio of tensor to scalar 
perturbations. As the potential for moving beyond first order perturbations has been explored, these three observable 
quantities have been joined by a measure of the departure from gaussianity exhibited by the perturbations, the non- 
gaussianity parameter /nl- This parameter is not yet well constrained by observational data in comparison with the 
other quantities but can already be used to rule out models with particularly strong non-gaussian signatures. 

There are two main approaches to studying higher order effects and non-gaussianity. One approach uses nonlinear 
theory and a gradient expansion in various guises, either explicitly, e.g. Refs. [3,4 or in the form of the AiV formalism, 
e.g. Refs. [3 EJ 13 HE EH E] By virtue of having to employ a gradient expansion this approach is so far only usable 
on scales much larger than the particle horizon. The other approach uses cosmological perturbation theory following 
Bardeen [H] and extending it to second order, e.g. Refs. [T3[II1[T33[T1[TJ[TE1^ (for 
an extensive list of references and a recent review on these issues see Ref. [28]). This approach works on all scales, 
but can be more complex than in particular the AiV formalism. Both these approaches give the same results on large 
scales [29 . We will follow the Bardeen approach in this paper. 

As the first order perturbations of the inflaton field are taken in the standard treatment to be purely gaussian it 
is in general necessary to go to second order in order to understand and estimate the non-gaussian contribution of 
any inflationary model (for a recent review see Ref. [28]). Deriving the equations of motion is not trivial at second 
order and only recently was the Klein-Gordon equation for scalar fields derived in closed form, taking into account 
metric backreaction |30j . This allows for the first time a direct computation of the second order perturbation in full, 
in contrast with previous attempts which have focused only on certain terms in the expression, for example Ref. |31j . 
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In this paper we solve numerically the second order Klein-Gordon equation in closed form in Fourier space and show 
that this procedure is readily applicable to the study of non-gaussianity and other higher order effects. As this is, to 
our knowledge, the first numerical solution to the full second order evolution equation we will outline the numerical 
steps taken in the system we have developed, examine the current constraints on the calculation and describe the 
next steps required in detail. This calculation uses the slow roll version of the second order equation, but solves the 
full non-slow roll equations for the background and first order. The models that we test in this paper are single field 
models with a canonical action. Significant second order corrections are expected only when a non-canonical action or 
multiple fields are used, or slow roll is violated. Numerical simulations will be particularly useful in analysing models 
with these characteristics. We will discuss in Section [V] planned future work to extend our current numerical system 
to deal with these extensions beyond the standard single field slow roll inflation. 

In Section [IT] we will give a brief outline of perturbation theory and describe the second order perturbation equations 
that will be numerically calculated. Section |III| describes the numerical implementation of the calculation, including 
the initial conditions used and the computational requirements. We present the results of this calculation in Section 
IV including a comparison of the second order perturbation calculated for the \m 2 (p 2 and |Ay? 4 potentials. We will 
cllscuss these results and the next stages of this work in Section |Vf 



Throughout this paper we set h = c = 1 and use the reduced Planck mass Mpl = v / 87rG. Overdots and primes 
denote differentiation with respect to our time variable n (the number of e-foldings) and conformal time rj, respectively, 
and will be defined explicitly when first used. We will work in a flat Friedmann-Robertson- Walker (FPW) background. 



II. PERTURBATIONS 



In this section we will briefly review the derivation of first and second order perturbations in the uniform curvature 
gauge and describe the slow roll approximation that we will use in this paper. There are many reviews on the subject 
of cosmological perturbation theory, and here we will follow Ref. [28 . The full closed Klein-Gordon equation for 
second order perturbations was recently given by one of the authors and we will outline the derivation in Ref. |30j 
below. 



A. First and Second Order 



In this paper we will consider perturbations of a single scalar field and will work throughout in the uniform curvature 
or flat gauge. Our goal is to describe scalar perturbations up to second order and the first step to achieve this is to 
examine the metric tensor: 

.g 00 = -a 2 (1 + 20! + 2 ) , (2.1) 

2 ( ~ 1 



90i - a 2 (Si + -B 2 j , (2.2) 
9ij = a 2 [(1 - 2Vi - ifa) Sa + 2E ltij + E 2Aj ) , (2.3) 

where a = a{rf) is the scale factor, rj conformal time, Sij is the flat background metric, 4>\ and (\> 2 the lapse functions, 
and "01 and ip 2 the curvature perturbations at first and second order; B\ and B 2 and E\ and E 2 are scalar perturbations 
describing the shear. Spatial 3-hypersurfaces are flat in our chosen gauge and so 

fa = fa = E x = E 2 = , (2.4) 

where the tilde denotes quantities in flat gauge. 

The Sasaki-Mukhanov variable, i.e. the field perturbation on uniform curvature hypersurfaces [32, 33 , evaluated at 
first order is given by 

SVi = S<Pi+ , (2.5) 

where (pa is the background value of the field and the perturbations of (p are defined as 

ip(x^) = <p (r]) + Spxiv, x l ) + ~S(p 2 {7], x l ) . (2.6) 
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At second order the Sasaki-Mukhanov variable becomes more complicated 

2 



c — c- fa i 
dip 2 = dip 2 + — ip2 + 
n 



2H<p' + <(% - ^V 



+ 2^Vi^i + - 2Si Pl , k E 1 k + X ty, E) 



n 



H 



(2.7) 



where X (tp, E) contains terms quadratic in gradients of the metric perturbations i/'i and E-±. From now on we will 
drop the tildes and talk only about variables in the flat gauge. The potential of the scalar field is also split 



U{tp) = U + 6Ui + ^SU 2 , 8Ux = U :V S Vl 



SU 2 = U jtp9 5tpi + U^6ip 2 , 



(2.8) 



where U m = 9U 
have 



ip — -q^. The Klein-Gordon equation describes the evolution of the scalar field. For the background field we 

(2.9) 



where 7i = — is related to the Hubble parameter H by H = aH. The first order equation is 

<W' + ZHStpi + 2a 2 U, ip <f> 1 - V 2 (5(pi - y' V 2 B x - ^ + a 2 U \ vlp 6<pi = , 

and the second order 

6ip 2 " + 2H5(p 2 ' - W 2 S(p 2 + a 2 XJ^8ip 2 + a 2 U \ vw {8ip{) 2 + 2a 2 U tV (j) 2 - cp' (W 2 B 2 + 4>' 2 ) 
+ VoBi,*^* + 2 (2Hip' a + a 2 U. v ) B lfk B^ + 40! (a 2 U,^5 Vl - V 2 ^i) + 4<p <f>i<t>i 



(2.10) 



2S<pi (V 2 Si 



(2.11) 



where as mentioned before all the variables are now in the flat gauge. 

The Einstein field equations are also required at first and second order. We will not reproduce them here but instead 
refer the interested reader to Section II B of Ref. [50]. Using the perturbed Einstein equations, the Klein-Gordon 
equations above can be written in closed form at both first and second orders. These equations will form the basis of 
the numerical scheme described in Section IIIII 

The dynamics of the scalar field becomes clearer in Fourier space but terms in the second order equation of the 



form (Sipi(x)) require the use of the convolution theorem (see for example Ref. 
we will write 5tp(k l ) for the Fourier component of 5tp(x) such that 

1 



(2tt) s 



d k5ip(k l ) exp(ifcjX J ) 



]). Following Refs. [3DJ and 

(2.12) 



where k l is the comoving wavenumber. 

In Fourier space the closed form of the first order Klein-Gordon equation transforms into 



Stptitf)" + 2H5tp 1 (k i )' + fc 2 5^i(F) + a 2 



8ttG 



2y' Q U >lp + (v'o) 



,8ttG, 
~H 



Stpxik*) = 0. 



(2.13) 



As mentioned above the second order equation requires more careful consideration with terms quadratic in the first 
order perturbation, which require convolutions of the form 

1 



f{x)g{x) 



(2tt) 2 



d 6 qd d p5 6 (k* - p l - q^WMq 1 ) . 



(2.14) 



For convenience we will group the terms with gradients of 6<pi(x) together and denote them by F. The full closed 
form second order Klein-Gordon equation in Fourier Space is 



+ 2nSip 2 '(k l ) + k 2 Sip 2 (k l ) + a 2 



8ttG 



2tp' U, 



,s 2 87tG 



SMk 1 ) 



(27T) 



d 3 qd 3 pS 3 (k l -p l -q l ) 



16ttG 

n 



[JO<^ V)<M<f ) + ^ a a ^w*vi(p i ) <y vi(g i )] 



/8ttG\ 



[■Jf) ^ [2a 2 U^' 5Mp l )$M<l l ) + ^X5^)5 Vl {q 1 )] 



2 f4nG\ 2 ^X 

V n ) h 

47rG" f _ f / i\ c ft i\ 2 

^^r~<fo d( pi ip )<Vi {q) + a 

-F(5ipi(k i ),S(f 1 '(k i )) = 0. 



IF 



8ttG 



S(p 1 (p i )5ip 1 (q l ) 



(2.15) 
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Here we use X — a 2 (8irGUo(p' /Tt + U tV ) for convenience. The F term contains gradients of Sipi in real space and 
therefore the convolution integrals include additional factors of k and q. It is given by 

F (S(f 1 (k'-),S(p 1 '(k i )) = JL J d 3 pd 3 q6 3 (k> -p* - A 2 ^fVi&f) (*<^l(<f ) + </o<V(<f )) 



p l qjk J k t 
k 2 



+2- (^Gy Pl q l Pm q m +p 2 q 2 
k 2 q 2 



J H\HJ 

4irG 



(p' 6<px(p l ) (XSip!(q l ) + tp'oStpiiq 1 )) 



H 



k 2 



V H H 



piq l p m g m 



H 



8nG 



k 2 



g <fyi(p )<Vi(<? ) - 



A- 2 



(2.16) 



B. Slow Roll approximation 



In order to establish the viability of a numerical calculation of the Klein-Gordon equation we have confined ourselves 
in this paper to studying the evolution in the slow roll approximation. In our case this involves taking 



(2.17) 



such that X — and H 2 — (8nG/3)a 2 Uo. The slow roll parameter e# as defined in Refs. [5U] and [25] (which is the 
square-root of the usual e) is given by 



H 



(2.18) 



With this approximation the second order equation (2.151 simplifies dramatically and with the F term included is 
dip 2 "(k l ) + 2H5<p 2 '(k l ) + k 2 5ip 2 (k l ) + (a 2 U w - 24^G(<^,) 2 ) 5<p 2 (k l ) (2.19) 



H- / <f'<l V - q'){'>- [ U tVVV + ^-tf/ Q U )VV j <Vi(p'><Vj(</) - 



16ttG 



8ttG 



d A p d A q5 6 {k l -p l -q l ) 



\n q 



-</>o 



, / / piq 1 +p 2 2 W'?' 



V fc 2 



1 q 2 +piq l 



k 2 



The numerical simulation in this paper will solve the slow roll version of the second order above, Eq. (2.19 1, the first 
order equation (2.131 and the background equation (2.9 1. In the next section we set up the correct form of these 
equations for the numerical simulation and discuss the implementation and some tests of the accuracy of the method. 



III. NUMERICS 



Our goal in this paper is to show that, just as at first order, a direct numerical calculation of the second order 
perturbations of a scalar field system is achievable and in this section we will outline how we have implemented this 
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system. In structuring the numerical system we have closely followed the work done at first order by Martin and 
Ringeval [36l[37] and previously by Salopek et al. 38J. 

A finite numerical range of k modes to be calculated is required. The upper cutoff in k, which marks the smallest 
scale considered, is well motivated by the difficulty in observing primordial perturbations at these small scales. At 
the other end we need to specify the largest scale or smallest k that we will consider. Analytically this is often taken 
to be the size of the universe, with k — being the equivalent mode. One immediate problem with this is that the 
Bunch-Davies vacuum initial conditions outlined in Section [ill B | blow up. The standard workaround is to implement 
a cutoff at large scales beyond which the amplitude of perturbations is zero. This is a pragmatic approach but recently 
there has been some evidence that a sharp cutoff similar to this could be responsible for the lack of power at large 
scales in the WMAP data [39l HOI HD H2] • 

The main concern is that the k range covers most if not all the modes observed to date in the CMB. The WMAP 
team rely for their main results, [2 , on £ multipoles in the range £ £ [3, 1000] which corresponds approximately 1 to 
k £ [0.92xl0" 60 ,3.1 x 10~ 58 ] M PL = [3.5x 10~ 4 , 0.12] Afpc" 1 . We will consider a similar range of k modes in this 
paper, taking three different ranges outlined in Section |TV| The choice of k range is flexible with the only constraint 
being that the number of modes at second order is one greater than a power of two. This enables faster integration 
using the Romberg method as explained below. 



A. Equations 



The equations in Section II B are not set up for a numerical calculation and in this section we rearrange them into a 
more suitable form. This involves a change of time coordinate and grouping of terms into smaller units for calculation. 
The second order slow roll equation (2.19) can be further simplified by performing the p integral and changing to 
spherical polar coordinates q,9,to where q = |q|. The d 3 q integral becomes 



<Tq 



q 2 dq 



2tt 



(3.1) 



For each k mode equation we take the 9 = 0,oj = axis in the direction of k l , so that the angle between k % and 
q l is 9 and the scalar product qik 1 — qk cos 9. The argument of each Sipi or dpi term depends on 9 through 
\k l - q l \ = y/k 2 + q 2 - 2kq cos 9 and so must remain inside the 9 integral. There is no uj dependence in Sipi with this 
choice of axes, so the last integral is simply evaluated. 

In the slow roll case there are only four different 9 dependent terms, here labelled A-D: 



A(fc\<f) = 


f sm(6)S(p 1 (k i - q^dd , 

Jo 




5(fc\«f) = 


/ cos(9)sm(9)6<p 1 (k l - q^dO , 
Jo 




C(fc\<f) = 


f sm(9)Sip 1 , (k t -q l )d9, 
Jo 




D(fcW) = 


[ cos(0) sin(0)<fy>i' \k l -q l )d9. 
Jo 


(3.2) 



Written using the terms in Eqs. (3.2 1 the slow roll equation (2.19) becomes: 

S<p 2 "(k l ) + 2H5^{k l ) + k 2 Si P2 (k t ) + (a 2 U^ v - 24^G(^) 2 ) 5ip 2 (k l ) + S(k l ) = , 
1 



(3.3) 



S(k l 



(27T) 



dq { a 2 U >ip(p(p q 2 6ipi(q i )A(k i ,q i ) 



8ttG 



H 



3a 2 U w + 7 -q 4 + 2k 2 q')A(k\q 1 )- { - 



k 2 



kq 3 B(k l ,q i ) 



8wG . 



-q 2 C{k\q l 



J,2 



kqDik^q 1 ) 



<W(9*) , 



(3.4) 



The approximate conversion for t is I ~ %- and a Megaparsec is given in Planck units as IMpc 1 ~ 2.6247x10 57 Mpl. 
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where S(k l ) is the source term which will be determined before the second order system is run. The full set of 
equations which must be evolved are then Eq. (2.9 1 for the background, Eq. (2.131 for the first order perturbations 



and Eqs. (3.3 1 and (3.4 1 for the second order and source terms. 



A more appropriate time variable for the numerical simulation is the number of e-foldings, and hence we use 

n = log(a/a in it) , (3.5) 



as our time variable instead of conformal time. Here, a in ; t is the value of a at the beginning of the simulation. If 
a is set to be 1 today we can calculate ai n j t once the background run is complete and the end time of inflation is 



determined as in Section III C We will use an overdot to denote differentiation with respect to n. 
The changes in derivatives required are as follows: 

d dn d d 

dn dn dn dn ' 

d dn dn d d 

dt dt dn dn dn ' 



(3.6) 
(3.7) 



where n and t are conformal and coordinate time respectively with H = % ja and Ti — aH. As mentioned above the 
value of a at the end of inflation is calculated using the connection equation (see for example Eq. (3.19) in Ref. [35] or 
Eq. (7) in Ref. |43p assuming that instantaneous reheating occurs at the end of inflation. This gives approximately 
65 e-foldings from the end of inflation until now. The background and first order equations written in terms of the 
new time variable n are 



U 



H 2 



0. 



H 
H 



6<pi 



k 



8ttG 



H 2 H 2 



6ip!=0. 



(3.8) 
(3.9) 



The second order equation in terms of n is 



<wm (3+|W 2 (fc*)+ (A 



8<P2 (fc l 



U 



ipip 



H 2 



247tG(<a,) 2 5ip 2 {k l ) + S(k l ) = . 



(3.10) 



S(k>) 



1 



(2n) 2 
8ttG 



dq 



(aH) 2 
- 8irG(po 



U 



j^ q 2 5 Vl {q l )A{k\ q >) 
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3a z U vv q' + -g 4 + 2k z q z A{k\ q l ) - - + ^ kq A B{k\q l ) 



9 q 2 



q 2 C(k\q*)+[2 l2 



l 2 )kqD(k 



1 q l ) 



(3.11) 



where 



Cik^q 1 ) = — Cltf-q 1 )^ / sm{6)5ipi(k l - q l )d6 . 
aH Jo 



D{k\q l ) = —D{k l -q* 



cos{6) sin(6)Sipi(k i - q^dO . 



(3.12) 



The argument of Sipi and Sipi in the A—D terms requires special consideration. To compute the integrals, 9 is 
sampled at 

N 6 = 2 l + l (3.13) 
points in the range [0, ir] (for some / G N to allow Romberg integration) and the magnitude of k l — q l is found using 



\k l - q l \ = y/k 2 + q 2 -2kqcos(6) 



(3.14) 
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While 5ipi(k l ) = 5(pi(k), the value of \k l — q l \ is at most 2fc max where k, q S [k m i n , A; max ]. This means that to calculate 
the source term for the k range described we require that 8ip\ and 8ip\ be known in the range [0,2£; max ]. In Section 



IIIC we will show that this first order upper bound does not significantly affect performance. On the other hand 
\k l — q l \ can also drop below the lower cutoff of calculated k modes. As discussed above we will implement a sharp 
cut off and take Sipi(k) = for the values below fc min . When Ak ~ fc m j n this affects only the k = q modes and is only 
significant close to fc m i n . Section [ill D | describes how the accuracy is affected by changing Ak and other parameters. 
Without extrapolating outside our computed k range it appears to be very difficult to avoid taking this small number 
of 5ifi s to be zero. 

The value of \k l — q l \ will not in general coincide with the computed k values of Sipi. We use linear interpolation 
between the closest fcs to estimate the value of S(fi at these points. We leave to future work the implementation of a 
more accurate but also numerically intensive interpolation scheme. 

Throughout the discussion above we have not specified any particular potential U and indeed the numerical code 
can use any reasonable potential provided that it gives a period of inflationary expansion in the e- folding range being 
simulated. In this paper we have used the two standard potentials U — ^m 2 ip 2 and U = \Xf 4 but a modular system 
allows another potential to be used instead. We choose the parameters to and A in agreement with the first order 
perturbation results from WMAP5 at the pivot scale fcwMAP = 0.002Afpc~ 1 ~ 5.25 x 1CP 60 Mpl with the values 
to = 6.3267 x lCT 6 Afp L , A = 1.5506 x 1CT 13 . 

B. Initial Conditions 

The background system requires initial conditions for ipo,<fo and H. These initial conditions and the range of 
c-foldings to be simulated must be selected with the choice of potential in mind. Not only must the e-folding range 
include an inflationary period, but the k modes to be calculated at first and second order must begin inside the 
horizon during this range. For example the initial conditions (pa — 18AfpL, </3o = — 1-WpLj H — 4.65x10 -5 Mpl for the 
\rv?^p 1 model give the background evolution described below and shown in Figure [lj 

The initial conditions are set for each k mode a few e-foldings before horizon crossing. This follows the example of 
Salopek et al. [38] and is justified on the basis that the mode is sufficiently inside the horizon for the Minkowski limit 
to be taken. This initial time, riinit(fc), is calculated to be when 

—^— = 50. (3.15) 
an | in it 

The range of e-foldings being used must include the starting point for all k modes, but the parameter on the right 
hand side, here chosen to be 50, can be changed if needed. We use the small wavelength solution of the first order 
equations as the initial conditions [38J, with 

V8nG e~ ikr> 

<fyl|init = /=-, (3.16) 

a y/2k 

, \/8jtGc^ A . k \ 



where the conformal time ij can be calculated from r\ = J dn/aH ~ —(aH(l — en)) 1 , when en changes slowly. For 
example fcwMAP is initialised about 65 e-foldings before the end of inflation and crosses the horizon about 5 e-foldings 



later. We also use these formulae in the calculation of the source term in Eq. (3.11 ) to determine the value of Scpi for 
a k mode before its evolution starts. 

We are interested in the production of second order effects by the evolution of the the gaussian first order modes and 
we make no assumptions about the existence of second order perturbations before the simulation begins. Therefore 
we set the initial condition for each second order perturbation mode to be 5<f2 = 0, 5<f2 = at the time when the 
corresponding first order perturbation is initialised. 

C. Implementation 

The current implementation of the code is mainly in Python and uses the Numerical and Scientific Python modules 
for their strong compiled array support [44] . The core of the model computation is a customised Runge-Kutta 4th 
order method (see for example Eq. (25.5.10) in [15]). Following Refs. [351 [37] the numerical calculation proceeds in 



four stages. The background equation (3.8), rewritten as two first order (in the time derivative) equations, is evolved 
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FIG. 1: The end of inflation is determined by calculating when eh = —H/H = 1 (red dashed line). Along the a:-axis, n is the 
number of e-foldings from the start of the simulation. 



from the specified initial state until some end time required to be after the end of the inflationary regime. The end 
of inflation occurs when d 2 a/dt 2 is no longer positive and the parameter £jj = —H/H first becomes greater than or 
equal to 1 (see Figure [TJ. Here, this specifies a new end time for the 1st order run, although the simulation can run 
beyond the strict end of inflation if required. The initial conditions for the first order system are then set as outlined 
above. 



The system of ordinary differential equations for the first order perturbations from Eq. (3.9 1 is calculated using a 



standard Runge-Kutta method. A fixed time step method is used in order to simplify the construction of the second 
order source term and because a "priori it is not known which time steps would be required at second order if an 
adaptive time step system were used. The first order equations are separable in terms of k and so it is straightforward 
to run multiple instances of the system and collate the results at the end. However, as will be discussed below, the 
first order calculation is not computationally expensive in comparison with the other stages and takes of the order of 
a few minutes for around 8000 time steps and 1025 k modes. 

Once the first order system has been solved the source term for the second order system must be calculated. As the 
real space equation for the source involves terms quadratic in the first order perturbation it is necessary to perform 
a convolution in Fourier space, as shown in Eq. (3.3 1. Transforming back into real space was not considered due to 
the presence of both gradient operators and their inverses. Here the slow roll version of the source term integrand 
has been used, but the method can equally be applied to the full equation. This stage is the most computationally 
intensive and can be run in parallel as each time step is independent of the others. The nature of the convolution 
integral and the dependence of the first order perturbation on the absolute value of its arguments requires that twice 
as many k modes are calculated at first order than are desired at second order as explained above. As the first order 
calculation is computationally cheaper than the source term integration, this does not significantly lower the possible 
resolution in /c-space, which is still limited by the source term computation time. Once the integrand is determined 
it is fed into a Romberg integration scheme. As for 6 which was discretised by Ng points in Eq. (3.2), this requires 
that the number of k modes is 



(3.18) 



for some 2 ^ £ N. This requirement can be lifted by opting for a less accurate and somewhat slower standard quadrature 



2 The number of discretised k modes JVj. does not need to be equal to Ng. 
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routine. 

The second order system is finally run with the source term and other necessary data being read as required from 
the memory or disk. The Runge-Kutta method calculates half time steps for each required point, for example if y(x n ) 
is known and y(x n +i) = y(x n + h) is required (for step size h), the method will calculate the derivatives of y at 
y{x n ),y(x n + h/2) and y(x n + h). As we need to specify the source term at every calculated timestep, the requested 
timcstep for the second order method must be twice that used at first order. This decreases the accuracy of the 
method but does not require the use of splines and interpolation techniques to determine background and 1st order 
variables between time steps. 

The second order system is similar in run time to the first order system but the source integration is more complex, 
involving the integration of N 2 x Ng values at each time step. Although a large amount of data is produced at each 
step at this stage for each of the wavenumbers fc, only the integrated result is stored to be used in the second order 
run. Results for each stage are stored in the open HDF5 standard which can deal efficiently with large files and is 
very portable, allowing data analysis independent of the Python/Numpy programming environment. We intend to 
release the program under a suitable license once the code has matured and some of the improvements discussed in 
Section [V] have been implemented. 



D. Code Tests 



We have tested the numerical code in a variety of controlled circumstances in order to quantify the effect of different 
choices of parameters. In particular it is important to know whether the values picked for Ng, the number of discretised 
8s, Afc, the size of the spacing of the discretised fc modes, and the range of fc values significantly impacts on the results. 
The ODE solving parts of the code are straightforward and follow standard algorithms. 

As mentioned above the WMAP results [5] use observations in the range k £ [0.92x 10 - 60 , 3.1 x 10~ 58 ]Mpl = 
[3.5xlCP 4 , 0.12] M^. 1 . We will consider three different k ranges both in our results and the tests of the code 3 : 



K Y = [l.9xl0~ 5 , 0.039] Mpc" 1 , Afc = 3.8xlO~ 5 Mpc 
K 2 = [5.71 x 10" 5 , 0.12] Mpc^ 1 , Afc = 1.2x 10~ 4 Mpc 
K 3 = [0.52xl0~ 5 ,0.39] Mpc" 1 , Afc = 3.8xlO" 4 Mpc 



(3.19) 



The first, K\, has a very fine resolution but covers only a small portion of the WMAP range. The next, K2, is 
closest to the WMAP range with a still quite fine resolution. The final range, K 3 , has a larger fc mode step size 
Afc = lxlO _60 MpL = 3.8 x 10 _4 Mpc _1 and covers a greater range than the others, extending to much smaller scales 
than WMAP can observe. 

The main new addition in the code is the calculation of the convolution of the perturbations for the source term 
Eq. (3.111. In particular the first of the dependent terms in Eq. (3.2), A, can be convolved analytically for certain 



smooth 5<pi(k)s. We take Sifi(k) to be similar in form to the initial conditions (3.161, for example 5<pi(k) cx 1/vfc 
with proportionality constant a. If I a denotes the convolution of the A term: 



I A {k) = 2ir 



q 2 6ip 1 (q)A(k,q)dq, 



(3.20) 



then putting in S(fi(k) = a/vk gives 



I A (k) = 2-ko? 



dqq^ d0(k 2 + q 2 - 2fc(j cos (9) " 4 sin 6» . 
Jo 



(3.21) 



3 The k ranges in AfpL are: 

K\ = [0.5xl0 -61 , 1.0245X10" 58 ] Af PL , Afc = lx 10" 61 M PL 

K 2 = [1.5 x 10" 61 ,3.0735xl0" 58 ] M PL , Afc = 3x 10" 61 M PL 

K 3 = [0.25xl0" 60 ,1.02425xl0~ 57 ] M PL , Afc = lxlO" 60 M PL . 
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(a)The relative error for different Ng, the number of discretised (b)The relative error for the 3 different k ranges KX, K2, K3 

8s, keeping the other parameters fixed and using the K3 range. (starting from the left). The parameter AA; is set equal to 

The upper blue line (Ng = 129) and middle green line lxlO~ 61 M PL , 3x 1CT 61 A/ PL , lxlCr 60 Af PL respectively. 
(Ng = 257) have relative errors at least an order of magnitude 
larger than the lower red line (Ng = 513). 

FIG. 2: Comparison of relative errors for different Ng and k ranges. 



This has the analytic solution 



A(k) 



18 fc 



3k 6 



log (2V*) 



l] 



7T 

2 



arctan 



+ log (2 (y/k~ + y/k 

min 1 



-fc)) - log (2(7^+^ n 

+ \Amax sj fc max - k (-3k 2 + 14fcfc max - 8fc max ) + l/fcmax 



k j 
It (3k 2 



'" 2 \ 



14fcfc max + 8k 2 

— V^^min \/ k — fc m in (3fc — 14fcfc max + 8fc, nax ) — y/ k + fc m j n (3k + 14fcfc max + 8fc max ) I . (3.22) 
Wc have tested our code against this analytic solution for various combinations of k ranges and Ng. The relative error 



I analytic — calculated! 
I analytic | 



(3.23) 



is small for all the tested cases but certain combinations of parameters turned out to be better than others. The 
relative error of all the following results is not affected by the choice of a so we will keep it constant throughout. 

We first tested the effect of changing Ng, the number of samples of the 6 range [0,ir). Figure 2(a) plots these 
results for the k range K3 with Afc = 1x10 _60 Mpl. Only three values of Ng are shown for clarity. It can be seen 
that increasing Ng decreases the relative error (for the convolution term at least) when the other parameters are kept 
constant, as one might expect. 

As mentioned above the choice of k range is especially important as the convolution of the terms depends 
heavily on the minimum and maximum values of this range. Indeed this is clear from the analytic solution in 
Eq. (3.221. Figure 2(b) shows the difference in relative error for the three different k ranges described above with 
3.8xl0~ 5 ,1.2xl0- 4 and 3.8x 10 _4 A/pc _1 (Afc = lxl0~ 61 , 3xl0~ 61 , lxl0" 6O M PL ) respectively. The accuracy 



Afc 



is similar in all three cases. 

Another important check is whether the resolution of the fc range is fine enough. Varying Afc can not be done in 
isolation of course, if the constraint for iV^, Eq. (3.18 1, is to be met. For this test the end of the fc range changed with 
Afc but the other parameters were kept fixed as fc min = lxlO~ 60 Af PL = 3.8x 10~ 4 Mpc~\ N k = 1025 and Ng = 513. 
Figure [3] plots these results again for only three indicative values. For Afc < fc m ; n , here the upper two lines, there is 
a marked degradation in the accuracy of the method. This is understandable as many interpolations of multiples of 
Afc below fc m i n will be set to 0. Once Afc is greater than fc m i n the relative error is very similar for higher values (not 
shown in the figure). 

It should be noted that these tests show the relative errors in the computation of the A convolution term, the most 



straightforward term in Eq. (3.2 1, only and do not represent errors for the full calculation. However, they show that 
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at least for the pure convolution term the accuracy is good compared with the analytic results. Equation (3.22) gives 
some indication of the difficulty involved in finding an analytic solution for the other terms, although this is a goal for 
future work. Having described the implementation and accuracy of the numerical system we will outline our results 
in the next section. 



IV. RESULTS 



The main result of this paper is the demonstration of a numerical solution to the closed Klein-Gordon equation of 
motion for second order scalar field perturbations as described in Eq. (2.19). This includes the slow roll approximation 



of the source term for second order perturbations, but we have not used a slow roll version of the evolution equations 
for the background or first order perturbations. 

As a proof of concept we have tested the system with two standard potentials, ^m 2 ip 2 and \\<p 4 and computed 
results across three different k ranges. As expected, considering the use of a single slowly rolling field, the second 
order perturbation we have calculated is extremely small in comparison with the first order term. However there 
are already differences apparent between the two potentials which will be outlined below. We have calibrated the 
parameters m and A of the potentials using the WMAP 5 normalisation at /cwmap = 0.002Afpc~ 1 = 5.25x 10~ 60 Af] 
[2]. We have outlined in Eq. (3.191 the three k ranges that we will use, 



PL 



K l = [l.9xl(T 5 , 0.039] Mpc -1 , Ak = 3.8xl0- 5 Mpc- 1 
K 2 = [5.71xl0~ 5 ,0.12] Mpc" 1 , Ak = 1.2x lO^A/pc" 1 
K 3 = [0.52x 10~ 5 , 0.39] Mpc -1 , Ak = S.SxlO^Afpc" 1 . 

Many of the results will be quoted for fcwMAP which lies in all three of these ranges. 

Given that the first order perturbations for the chosen potentials give an almost scale invariant power spectrum 
with no running, it is no surprise that the results from the three different k ranges are very similar. The second order 
source term is somewhat dependent on the lower bound of k (upper bound on size). This is expected and in the scale 
invariant case a log divergence can be shown to exist 39J. We have implemented an arbitrary sharp cutoff at A; m j n 
below which 6<p\ is taken to be zero. As mentioned above there is some evidence to suggest that a similar cutoff is 
supported by the WMAP data [5T1I52]. 

At first order our solutions agree with previous work [3HJ [371 [35], with oscillations being damped until horizon 
crossing (when k = aH) after which the curvature perturbation becomes conserved. Figure [5] shows the real and 
imaginary parts of the first order perturbations from when the initial conditions are set at k/aH = 50 to just after 
horizon crossing. The x-axis for most of the following figures shows the number of e-foldings left until the end of 
inflation instead of the internally used time variable n. 

In Figure [6] we show the evolution of the second order perturbations for wavenumber /cwmap- As mentioned above 
the overall amplitude of the second order perturbations is many orders of magnitude smaller than the first order ones. 
In Figures [H] and [H] the field values have been rescaled by k 3 ^ 2 /(t/2tt) to allow a better appreciation of the magnitude 
of the resulting power spectra. 

The source term S(k l ) is calculated at each time step using the results of t he fir st ord er and backgro und r uns. This 
term drives the production of second order perturbations as shown in Eqs. (2.191 and (3.101. Figure 7(a) shows the 
absolute magnitude of the source term for a single k mode, fcwMAP, for all time steps calculated. Figure 8(a) shows 
how the source term changes with the choice of k range. After horizon crossing the source terms are approximately 
equal. Before horizon crossing however there is a strict hierarchy with the smaller k ranges, K\ and K 2 , leading to 
smaller source contributions. As stated in Section [III D[ Ak should be at least as large as fc m i n in order to reduce the 
error to a minimum. 

The source term is large at early times, and closely follows the form of the spectrum of the first order perturbations 
as can be seen from Figure |7(b)| It is useful to compare the magnitude of the source term with the other terms in 
the second order evolution equation (3.101. If we let T denote the other terms, 



T(k i )= (3+f ) 



k 



H 2 



2AnG{^) 2 Sip 2 (k l ) 



(4.1) 



then Figure 9(a)| shows the absolute magnitude of both S and T. It is clear that the source term is of comparable 
magnitude only early in the simulation. Figure [9(b) | shows a comparison of |5|/|T| for three different k values. The 
larger the k mode the closer in amplitude S is to the rest of the terms in the ODE. A priori it is not known where 
S will be large for a particular chosen potential and mode but once determined it could be possible to significantly 
reduce the time required for the simulation by only calculating S in the regions where it is important. 
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FIG. 5: The first order perturbation Sipi rescaled by k 3 ^ 2 /(\/2ir) from the beginning of the simulation until around horizon 
crossing (red dot-dashed line). The real (blue) and imaginary (green dashed) perturbations are shown for /cwmap- 



All the results quoted so far are for the |m 2 y 2 
m 2 ip 2 . Figure 



10 



compares the models for fc 



WMAP 



model. We have also tested the jXip 4 model and compared it to 
shows how the source term for the \X(p 4 model 



Figure 10(b) 



is larger than the one for ^m 2 ip 2 to begin, but crosses over after a few e-foldings. After horizon crossing the \\<p A 
source term is again larger. As the results at first order for both models are so similar it is to be expected that the 
second order perturbations would be closely related. 

In Figure 11 the value of |5| at the start of the evolution of 8<f2 for each k mode is shown. The magnitude of the 
source term is much smaller for larger ks (smaller scales). Because the smaller fcs begin their evolution earlier the 
relative difference in |5| is not as pronounced when measured at a single timestep (see for example Figure [8(b)] ) . It 
should also be remembered that the magnitude of other terms in the second order ODE is small for larger ks as shown 
by the ratio |5|/|T| in Figure [9] where T is defined above in Eq. (4.1|. 

The source term for all fcs can also be compared for different timesteps. In Figure [12] the upper blue line shows 
IS^fc)! around 69 e-foldings before the end of inflation when 6y>2 has been initialised for only the very smallest k 
modes. The middle green line shows |5| when all 8<p2 modes have been started. Finally the lower red line plots |5| 
after all modes have exited the horizon, around 52 e-foldings before the end of inflation. 



V. DISCUSSION AND CONCLUSION 



In this paper we have described the numerical solution of the evolution equations for second order scalar perturba- 
tions, using the closed form of the Klein-Gordon equation, Eq. (2.19 1. We demonstrate that direct calculation of field 
perturbations beyond first order using perturbation theory is readily achievable, though not trivial. 



For this first demonstration we have limited ourselves to considering the slow roll source term in Eq. (2.19) but 
without imposing slow roll on the evolution terms of the ODEs. We have investigated two standard potentials, |m 2 </> 2 
and |A(/> 4 , to demonstrate the capabilities of the system. The singularity at fc = which arises as larger and larger 
scales are considered is avoided by implementing a cutoff at small wavenumbers below fc m i n . This is a pragmatic choice 
necessary for the calculation, but as mentioned above there is some evidence that such a cutoff might also explain 
lack of power at large scales in the WMAP data [101 EH H2] • It is also necessary to pick a maximum fc value, and this 
choice is dictated by computational resources and with reference to observationally relevant scales. In this paper we 
have used fc ranges which are comparable with the scales observed by WMAP. By comparing the analytical results 
of the convolution integral with the numerical calculation, we have chosen values of the parameters Ng,Ni~ and Afc 
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FIG. 6: The real (blue line) and imaginary (green dashed) components of the second order perturbation 8tp2 (fcwMAp) from the 
beginning of the simulation until around the time of horizon exit (red dot-dashed line). 
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FIG. 7: Source term and power spectrum for the WMAP pivot scale fcwMAP- 



which minimise the numerical error. The convolution scheme that we have implemented works best when Afc > fc m i n . 

We have shown explicitly that the second order calculations for our chosen potentials are obtainable once the cut-off 
for fc m i n is implemented. As expected for these unexceptional potentials in the slowly rolling regime the magnitude 
of second order perturbations is extremely suppressed in comparison with the first order amplitude. We have shown 
the evolution of the source term equation during the inflationary regime can be readily calculated. 

There are many possible next steps to improve the program outlined in Section |III| Chief amongst these is to 
implement the full source term equation (3.111. Although clearly more complicated than the slow roll case in Eq. (3.10 1 
only two more 9 dependent terms need to be added to A-D in Eq. (3.2 1. For the two test models we have used in 
this paper, which are both slowly rolling during inflation, it is not expected that using the full source equation would 
result in an appreciably different outcome until the end of the inflationary phase. Though once the field has stopped 
to roll slowly, new observable features might arise as is indeed the case at first order. 
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FIG. 8: Two different comparisons of the source term S. 
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(a)The source term (lower blue line) is compared with the T (b)The quotient of S and T terms for three different k values 
term (upper green line) as defined in the text for fcwMAP- The including the WMAP pivot scale. Depending on k the source 
source term is of comparable magnitude at the beginning of the term only dominates at early stages or is important throughout 
simulation. the evolution. 

FIG. 9: Source term 5* compared with T. 



Beyond this the next likely step is to implement a multi-field version of the system. This would allow the investi- 
gation of models that inherently produce large second order perturbations. In Ref. |30] the Klein- Gordon equation is 
given for multiple fields and upgrading the simulation to use these equations is a straight-forward if lengthy process. 

The performance of the numerical simulation could also be improved by analysing the most time consuming processes 
and investigating what optimisations could be implemented. As we have discussed above we have set = 1025 for our 
test runs. This provides good coverage of the WMAP k range but it is not clear whether it sufficiently approximates 
the integral to infinity for the source term. Currently we are restricted in our choice of iVfc by logistical factors i.e. 
the running time and memory usage of the code. By optimising the routines for both memory and speed it is hoped 
we can extend the maximum value of k to larger values. 

By computing the perturbations to second order we have direct access to the non-gaussianity of 5(p . While useful for 
the toy models discussed above (with /nl — 0) , when used to investigate models with predictions of large non-linearity 
parameter /nl this technique could yield greater insight into the formation and development of the non-gaussian 
contributions by studying the contribution of the different terms in the source term Eq. (3.11 ). It was shown recently 
that in order to calculate /nl instead of using the standard method based on the Lagrangian formalism [5U], one can 
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simulation, before horizon crossing. 



FIG. 10: A comparison of the source term for the \rri l tp l and j\<p 4 models. 




FIG. 11: The absolute magnitude of the source term at the initial start time for each k when k = aH x 50 deep inside the 
horizon. The results are for the range K x = [l.9x 10~ 5 , 0.039] Mpc" 1 = [0.5xl0" 8 \ 1.0245xl0~ 58 l M PL . 



instead use the field equations [351 HZ] • The method presented here will therefore eventually allow a full numerical 
calculation of /nl- 

In summary, we have demonstrated that numerically solving the closed Klein-Gordon equation for second order 
perturbations is possible. We have used the slow roll version of the source term in this paper, but hope to extend 
our work to use the full source soon. The two test models we have used have been shown to have negligible second 
order perturbations in line with analytic results. We have compared the analytic and numerical solutions for the 
convolution term and found them to be in good agreement. 
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FIG. 12: The absolute magnitude of the source term for all ks in the range at three different timesteps: the upper blue 
line when only the largest modes have been initialised; the middle green line when all modes have been initialised; and the 
lower red dashed line when all modes have exited the horizon. The k range shown here is Ki — [l.9x 10 _J , 0.039] Mpc -1 = 
[0.5xl0" 61 ,1.0245xl0~ 58 ] M PL . 
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